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Abstract 

In this work, we develop a probabilistic estimator for the voltage-to-current 
map arising in electrical impedance tomography. This novel so-called par¬ 
tially reflecting random walk on spheres estimator enables Monte Carlo meth¬ 
ods to compute the voltage-to-current map in an embarrassingly parallel 
manner, which is an important issue with regard to the corresponding in¬ 
verse problem. Our method uses the well-known random walk on spheres 
algorithm inside subdomains where the diffusion coefficient is constant and 
employs replacement techniques motivated by hnite difference discretization 
to deal with both mixed boundary conditions and interface transmission con¬ 
ditions. We analyze the global bias and the variance of the new estimator 
both theoretically and experimentally. In a second step, the variance is con¬ 
siderably reduced via a novel control variate conditional sampling technique. 

Keywords: Monte Carlo methods, electrical impedance tomography, 
random walk on spheres, reflecting Brownian motion, discontinuous 
diffusion coefficient, random diffusion coefficient, variance reduction 


1. Introduction 

The mathematical formulation of static electrical impedance tomography 
(EIT) leads to a nonlinear and ill-posed inverse problem, which is unstable 
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with respect to measurement and modeling errors. Namely the reconstruction 
of the real-valued conductivity n in the elliptic conductivity equation 

V ■ {kS/u) = 0 in D (1) 

from boundary measurements of the electric potential u and the correspond¬ 
ing current on the boundary of a bounded, convex domain D C d = 2,3, 
with piecewise smooth boundary dD and connected complement. Due to 
the limited capabilities of static EIT, many practical applications focus on 
the detection of conductivity anomalies in a known background conductivity 
rather than conductivity imaging, cf., e.g., Pursiainen [2H] and the recent 
work |12] by the second author. In this work, we consider such an anomaly 
detection problem, where a perfectly conducting inclusion occupies a region 
T inside the domain D. A possible practical application modeled by this 
setting is breast cancer detection, where the electric conductivity of high- 
water-content tissue, such as malignant tumors, is approximately one order 
of magnitude higher than the conductivity of low-water-content tissue, such 
as fat, which is the main component of healthy breast tissue, cf. [8]. 

The most accurate mathematical forward model for real-life impedance 
tomography is the complete electrode model (CEM), cf. [H], where the elec¬ 
tric potential u is assumed to satisfy the Robin boundary condition 

zp ■ KVu\dr) + u\qi:, = (f) ondD. (2) 

Here v denotes the outer unit normal vector on dD and the positive constant 
2 ; is the so-called contact impedance which accounts for electrochemical effects 
at the electrode-skin interface. Given the full Robin-to-Neumann map 

Rz,n : 0 H- 1/ ■ kVu\od, 

that maps the potential on the boundary to the corresponding current across 
the boundary, this knowledge uniquely determines 2 : and is hence equivalent 
to the knowledge of the Dirichlet-to-Neumann map. In this case, uniqueness 
of solutions to the inverse conductivity problem for isotropic conductivities 
has been proved under various assumptions on both, spatial dimension and 
regularity of the conductivity, cf., e.g., the works by Astala and Paivarinta 
|2] for d = 2 and Haberman and Tataru |T9] for d = 3. 

Notice that the operator ^ corresponds to idealized measurements on 
the whole boundary dD. In practice, however, only a hnite number of hnite- 
sized electrodes is available and thus only incomplete and noisy measurements 
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of the Robin-to-Neumann map can be obtained. Given snch discrete voltage- 
to-current maps, the nse of a regnlarization strategy is mandatory becanse 
of the severe ill-posedness of the inverse problem, cf. Alessandrini [1]. In 
statistical inversion theory, the inverse problem is therefore formnlated in 
the framework of Bayesian statistics, that is, all the variables inclnded in the 
mathematical model are treated as random variables. The solntion to the 
statistical inverse problem is then given by the posterior probability distri- 
bntion of the nnknown parameters conditioned on the measnred data, see 
e.g., Compnting the conditional mean estimate as well as com¬ 

mon spread estimates from the posterior density leads to high-dimensional 
integration problems and Markov chain Monte Carlo (MCMC) techniqnes 
are nsnally employed for this task. However, each sampling step in snch an 
algorithm reqnires solving the forward problem Q, (|^ nnmerically so that 
the compntation time can easily become excessive. This effect is amplihed by 
the fact that the Robin bonndary condition (|^ leads to singnlarities of the 
solntion u at the end points of the electrodes snch that nnmerical approxima¬ 
tions, both via hnite element and bonndary element methods, reqnire very 
hne discretization. 

In this work, we are concerned with the forward problem of EIT. More 
precisely, we develop a probabilistic estimator for the voltage-to-cnrrent map 
which has potential to overcome the aforementioned drawback if it is used 
on massively parallel hardware, such as GPUs, within the so-called Bayesian 
modeling error approach, cf. Kaipio and Somersalo 1201. The main advan¬ 
tage of the proposed method, beside its inherent parallel scalability, comes 
from the fact that the error estimates required for the Bayesian modeling 
error approach may be computed adaptively and on the fly at almost no 
additional computational cost. On top of that, our approach is well suited 
for uncertainty quantihcation in problems with random parameters. 

Due to the advent of multicore computing architectures, probabilistic es¬ 
timators for the numerical solution of boundary value problems for PDE in 
three or more dimensions have become a valuable alternative to deterministic 
methods. This is particularly true, when one needs to compute the solution 
at only a few points, or when moderate accuracy is sufficient. For instance 
in biophysical applications, where the linearized Poisson-Boltzmann equa¬ 
tion must be solved, efficient probabilistic numerical algorithms have been 
developed recently, see e.g. IMIIS]. However, in contrast to these works, the 
derivation of a probabilistic estimator for the voltage-to-current map corre¬ 
sponds to the approximation of paths of the partially reflecting Brownian 
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motion, cf. [TS], rather than the killed Brownian motion. The partially re¬ 
flecting Brownian motion behaves like the standard Brownian motion inside 
the domain and it is prevented from leaving the domain either by absorption 
or by instantaneons reflection. Under qnite general assnmptions, a Feynman- 
Kac type representation formnla in terms of the bonndary local time process 
of the partially reflecting Brownian motion for the electric potentials in FIT 
was recently obtained by Piiroinen and the second anthor in |37] . It is, how¬ 
ever, well-known that direct simnlation of the nnderlying Lebesgne-Stieltjes 
integrals with respect to the bonndary local time process is qnite a difficnlt 
task, see e.g. P HH [IS]- To be precise, the first order convergence obtained 
by Gobet’s half-space approximation scheme [IS] is cnrrently the state of the 
art. 

In this work we propose a different approach, namely we discretize with 
respect to space by expressing the nnknown electrical potential as the expec¬ 
tation of some anxiliary random variable obtained via a local finite difference 
discretization. This yields a novel second order space discretization scheme. 
A similar techniqne, nsing a first order approximation, was first introdnced 
by Mascagni and Simonov in [33] in the context of simnlation of diffnsion pro¬ 
cesses in discontinnous media. Also for the simnlation of diffnsion processes 
in discontinnous media, second order schemes were proposed and analyzed 
by Bossy et ah jS] and by Lejay and the first author [23]. The idea to use 
a local finite difference discretization for the simulation of the boundary be¬ 
havior of reflecting diffusion processes was introduced recently by the first 
author and Tanre [30] and further developed by the first author and Nguyen 
|31j . Other, related schemes were defined by Lejay and Pichot [26] and Lejay 
and the first author [2H]- The main advantage of the method proposed in 
this work, in comparison to the aforementioned works, lies in the fact that 
the variance of our estimator is greatly reduced due to a combined control 
variates conditional sampling technique. Therefore, we expect the method to 
become a valuable alternative to established deterministic methods for the 
problem at hand. 

The rest of the paper is structured as follows: We start in Section 2 by 
describing briefly the modeling of electrode measurements and the anomaly 
detection problem in FIT. In Section 3 we recall the basic idea of the random 
walk on spheres (RWOS). Subsequently in Section 4 we introduce the novel 
partially reflecting random walk on spheres estimator and in Section 5 we 
describe how it can be used to approximate electrode measurements. Then 
in Section 6 the variance reduction technique is explained. Section 7 gener- 
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alizes the partially reflecting random walk on spheres estimator to problems 
with layered conductivities as well as problems with random parameters. In 
Section 8 we present numerical examples to illustrate the efficiency of our 
algorithm. Finally, we conclude with a brief summary of our results and 
comment on directions for future research. 

2. Modeling of electrode measurements 

Let us briefly recall both, the modeling of electrode measurements in FIT 
and the anomaly detection problem. Consider the time-harmonic Maxwell’s 
equations, to be precise, Faraday’s law and Ampere’s law 

curlF = curliL = —{ioje — k)E, 

where u is the frequency, p the magnetic permeability and £ the electric 
permittivity. Notice that the physically relevant electric and magnetic held 
are given by the real parts and respectively. Now 

let us specialize on the static, respectively quasi-static case, i.e. direct input 
currents, respectively low frequencies u. Then the imaginary part of the 
electrical admittivity iue — k becomes negligible as well as the term iu^H. 
It can indeed be shown that the Maxwell system is approximated by 

curlF = 0, curlif = kE, (3) 

see e.g. [B]. In particular the electric held must be a gradient held E = — Vu 
for the scalar electric potential u. Substitution of this expression into the sec¬ 
ond equation in (|^ and taking the divergence Anally yields the conductivity 
equation Q. 

In anomaly detection problems, it is commonly assumed that the electric 
conductivity is constant apart from the anomaly. Without loss of generality 
let us assume that k = 1 in D\T. Moreover, we assume that T is simply 
connected and has a smooth boundary dT. In the setting which we are 
interested in here, T models a perfect conductor. Then potential differences 
in T equalize instantaneously and the governing conductivity equation is the 
Laplace equation 

Am = 0 in D\T (4) 

with Dirichlet boundary condition on dT 

u\dT = C, (5) 
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where the constant c is implicitly defined through the imposed electrode volt¬ 
ages. We consider here discrete voltage-to-current measurements performed 
using N electrodes Ei,...,En, attached to dD. The electrodes are modeled 
by disjoint surface patches which are assumed to be simply connected subsets 
of dD, each having a smooth boundary curve. Within the CEM, given N 
electrodes, the electric potential u satisfies the Robin boundary condition 

zu ■ Vu\oD + fuloD = g on dD, (6) 

where the functions /, g : dD —)■ M are given by 

N N 

fix) ■=^Xi{x), g{x) :=^UiXiix). 

1=1 1=1 

Here, Xii') denotes the indicator function of the l-th electrode Ei and the vec¬ 
tor U = {Ui,..., Un)'^ denotes a prescribed electrode voltage pattern. Notice 
that the CEM accounts for two important physical phenomena: First, the 
shunting effect of the highly conducting electrodes and second the fact that 
the current densities are limited by the contact impedance z : dD —M+, 
which is caused by a thin, highly resistive layer at the electrode-skin interface. 
We always assume that the ground voltage has been chosen such that 

N 

J2Ui = 0. ( 7 ) 

1=1 

For a given voltage pattern U G satisfying Q, the equations Q, ([^ 
and ([^ uniquely define the potential-current pair (m, J) G H^iD) with 

electrode currents 

^ ■^u\aDd(T{x), l = l,...,N, 

\^l\ Jei 

satisfying the conservation of charges condition 

N 

E = 0- («) 

1=1 

cf. [H]. For simplicity of the presentation let us assume that \Ei\ = \E\, 
I = 1, ...,N, throughout this work. 
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Figure 1: Current density kS/u (arrows), equipotential lines and prescribed electrode volt¬ 
ages on _E 4 and E-j. The box is a perfectly conducting inclusion in unit background con¬ 
ductivity. The forward problem of EIT is to determine the electrode currents (Ji,Jio)^- 

Figures and [^illustrate the EIT forward problem using the CEM. Both 
hgures are computed using synthetic measurement data simulated via a £- 
nite element discretization, cf. m The values of the contact impedances 
are comparable to those measured in real-life impedance tomography, cf. [7] . 
It has been shown experimentally that the CEM can predict EIT electrode 
measurements up to measurement precision, cf. [HI Si]- In Figure [^notice the 
peaks of the current density near the electrode edges caused by the shunting 
effect, which lead to severe difficulties in numerical approximations via de¬ 
terministic methods. In fact, the regularity of the potential decreases as the 
contact impedance tends to zero, cf. [10], which is a huge drawback since in 
practice one typically aims for good contacts, i.e., small contact impedances. 


3. The standard random walk on spheres 

The random walk on spheres (RWOS) estimator is a classical tool in 
stochastic numerics for elliptic and parabolic boundary value problems, orig¬ 
inally designed to solve the Dirichlet problem for the Laplace equation 

An = 0 in H, n = 0 on dD^ (9) 

see e.g. [33 E 2 i sni m]. For convenience of the reader and to introduce 
notations let us briefly recall the basic idea. Let xq E D and let doixo) 
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Figure 2: Boundary current density v ■ KVu|aD in the CEM corresponding to the setting 
of Figure The ticks on the x-axis correspond to the electrode midpoints. 

denote the radius of the largest sphere entirely contained in D and centered 
in xq. Then the classical Feynman-Kac representation formula, cf. [39], yields 

v(xo) = IE[n(lTr(5(xo,do(xo))))|Wd) = a;o], 

where W is the standard d-dimensional Brownian motion and t(S(xo, dnixo))) 
its hrst exit time from the sphere S(xo, dn(xo)), cf. [22]. As this representa¬ 
tion is valid for all points inside the sphere one may use the strong Markov 
property of the Brownian motion to obtain the conditional expectation 

'^(^o) = IE[n(fhV(5(a:i,d(a;i)))) |hfo = ^T{S{xo,dr>{xo))) = ^^i]. 

Due to the isotropy of the Brownian motion the points Wr{s{x,dD{ 3 :))) 
formly distributed over the sphere S'(x, ^^(x)), i.e. the above procedure de- 
hnes a time-homogenious Markov chain with state space {D, B{D)) 

and initial distribution given by the Dirac measure concentrated at Xq. 
The states of this chain can be computed via the relation 

Xj = Xj_i -|- Rjdoi^j-i), j > 1, 

where {Rj}j^fq is a sequence of random independent and isotropic unit vec¬ 
tors. For a particular realization of the random vector Xj we will use the 
lower case symbol Xj{u), where u is an elementary element from the probabil¬ 
ity space (D, R, Pa;^) of the Markov chain. We will suppress the to in our nota¬ 
tion if this causes no confusion. It can be shown that hmj_,.oo Xj = Xoo G dD, 














P 3 ;Q-a.s. and the value v{xoo) = 0(xoo) is given by the Dirichlet boundary 
condition, which yields the probabilistic estimator 

V : D X —)■ M, v{xo,Ljj) = v{Xoo{oj))- 

To obtain a practically realizable estimator, one usually introduces an 
£-layer 

De := {x G -D ; d{x, dD) < e}, 

where (i(-, dD) denotes the Euclidean distance to the boundary. Let So{xj, e) 
denote the surface of the sphere S{xj, doi^j)) that belongs to D^. The prob¬ 
ability of Xj+i lying in is then given by 

So{xj,e)T{d/2) 

27r^/^dD{xj))^-^' ^ ’ 

In the case of the Dirichlet problem we are interested in the discrete hrst 
hitting time, i.e., the index t{D^) = inf{j : Xj G D^}. By the spherical mean 
value theorem we have for all I G No 


E[n(X/+i)|xi, = E[n(X/+i)|xi] = n(xi), 

that is, the sequence {n(Xj)}jgNo is a discrete-time martingale with respect to 
{XjljgNg and thus by Doob’s optional stopping theorem, cf. [22], the stopped 
chain {n(xo),..., n(X.r(DE))} is a one as well, implying 

E[h(xo, ...,Xi,t{D^) >l]= E[v(Xr(D,))\xi,T{D^) >l]= v{xi). 

In particular we have E['0(xo, ■)] = and thus ^(xo, ■) is an unbiased 

estimator for the solution of the Dirichlet problem. The corresponding prac¬ 
tically realizable estimator is given by 


h^(xo,-) = 0(7ra^(X^(B^))), 

where 'T^dD{xr{De)) denotes the normal projection on the boundary dD. 

By the law of large numbers, n(xo) may be approximated by simulation 
of i.i.d. sample paths of the chain {Xi,..., Xt-(£)^)}. Moreover, one can show 
that the bias of the practically realizable estimator 'O'^(xo, •) is of order 0{e), 
£ —)■ 0, i.e., for sufficiently small £ there exist a constant C > 0 such that the 
root mean square error can be estimated by 


M 


E 

cf. |32lil]. 


l\M 


'^V%Xo,UJm) -V(xo) 


m=l 


/ ^^,2 , Var[r(xo,-)] 
+ M 
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4. The partially reflecting RWOS estimator 

Now let us turn to the derivation of the partially reflecting RWOS es¬ 
timator for the electrode currents Ji, I = For simplicity of the 

presentation we restrict ourselves here to the case d = 2 , however the gener¬ 
alization to d = 3 is straightforward. Throughout this section, we assume for 
simplicity of the analysis of the proposed estimator that the electrodes cover 
the whole boundary dD, more precisely, we consider the boundary condition 
(|^ with smooth function (j). The case of a finite number of discrete electrodes 
not covering the whole boundary is treated in the subsequent section. 

Let us consider the mixed Dirichlet-Robin boundary value problem arising 
from the anomaly detection problem Q, ([^, ([^ and let us assume for the 
moment that the constant c in ([^ is known. Moreover, we assume that 0 
is a smooth function. Let e be sufficiently small such that the e-layers do 
not intersect, i.e., = 0. In order to derive a probabilistic estimator 

for the potential u{x) at an arbitrary point x G D\T, we must approximate 
the partially reflecting Brownian motion starting in x with absorption in 
T. Therefore, let us define a time-homogeneous Markov chain 
with state space {D\T U {9}, Bq{D\T)), where we have adjoined an isolated 
cemetery point {9}. This cemetery point captures the missing mass, thus 
accounting for the fact that the chain is neither purely absorbing nor purely 
reflecting. The lifetime of the chain is (, that is, Xj = d for all J > C- 
By the strong Markov property of the Brownian motion we may use the 
standard RWOS, as long as the chain has not entered any of the e-layers. 
Now let xk denote an arbitrary state of the Markov chain inside one of the 
e-layers or T^, respectively. If xk G T^, the chain terminates and we have 
u{xk) = c -f- ro(a;x), where ro(x) = 0{e), e —)■ 0, for all x G T^. On the 
other hand, if xk G D^, the value of u{710d{xk)) is unknown. Without loss 
of generality let us assume that i^i'XdDixx)) = ei, where ei is the unit vector 
in direction (1,0)^. We consider a standard 5-point stencil finite difference 
approximation with stepsize h > 0 of the Laplacian 

A’^uIxk + hei) = ^ (u(xk) + u(xk + 2hei) + u(xk + hiei - 1 - 62 )) 

^ ( 11 ) 

u{xk + h(ei - 62 )) - ^u{xk + hei) j 

together with the second order one-sided finite difference approximation of 


10 


the normal derivative 


\/^u{xk) = “ ^u{xk) - u{xk + 2 hei ) j . 

If one of the points involved lies outside of D, we reduce h until all the points 
lie inside. Due to the boundary condition ([^ and the fact that u is harmonic 
in D\T we obtain thus 


u{xk) + zV'^u(xh) ^ + ri{xK)y 


( 12 ) 


where ri(x) = 0{h^ + e/h) for all x G Dg. Now we multiply equation by 
—2zh‘^ and equation (12) by 2h and sum them up which yields for the value 
u{xk) the following expression: 


uix^) = + hn{x^), 


h + z 


h + z 


( 13 ) 


where 


Rhu{xK) ■= ^ {u{xk + h(ei + Cs)) + u{xk + h(ei - Ca))). 


The key observation is that the expression (13) yields a probabilistic inter¬ 
pretation, namely the hrst term is the expected value of a random variable 
that takes the values u{xk + h(ei — ea)) and u{xk + h(ei -1- ea)) with equal 
reflection probability 


Pt{xk) = 


2{h -|- z) 

and the value 0 with absorption probability 

PaixK) = 


(14) 


(15) 


h + z 

We follow the approach pursued in [H] and recast this observation into 


an inhomogeneous integral equation of the second kind: 


u{x) = 
with the Radon measure 


D\T 


u{y)k{x,(iy) + F{x)^ x e D\T, 


{l-Pa{x))Ph(e^±e 2 ){x,dy), X E 

k{x, dy) := <( 0, x E 

P^\j,{x,dy), else, 


( 16 ) 
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and inhomogeneity 


s(x) + hri{x), 

F{x) := { u{7raT{x)) +ro(x), 

0 , 


xeD, 

X eT, 

else. 


/b Ph(ei±e 2 }i.x, dy), X G D\T, B G B{D\T) is the probability transition kernel 
corresponding to the random reflection in according to (13), whereas 


Pjj\^rp(x,dy), X G D\T, B G B{D\T) denotes the probability transition 
kernel corresponding to the standard RWOS in D\T. Finally, the score 
function of the random walk estimator is given by 


s(x) := XdA^) 


h(j){TlQD{x)) 
h + z 


( 17 ) 


In order to obtain a probabilistic estimator, we define a randomized ver¬ 
sion of the snccessive approximation of the partial snms of the Nenmann 


series 


FM 


F(x„), 

j=0 


where K denotes the integral operator in (16). A canonical choice for the 


probability transition kernel of the nnderlying Markov chain is 

obviously given by 

P{x,B) ■= [ k{x,dy), X eD\TU{d}, B eBa(D\T). (18) 
J B 

Let us denote the probability space of the Markov chain with transition 
kernel P given by (18) and initial distribution Xq ~ 5^0, Xq G D\T by 


(n, Pj;^) and let IExo['] denote the expectation with respect to the measure 

Pa;^. Then we may define the following partially reflecting RWOS estimator 
for the electrical potential 

u : D\T X n —)■ M 

as well as the corresponding practically realizable estimator 


u{xo,u) := y^^P{xj{u)) 

j=0 


(19) 


W : D\T X n —)■ M, u^{xo,u}) : = E s(a;j(a;)) + cxTflxt;-i{u)). (20) 


j=0 
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By the following result, both estimators (19) and (20) are convergent and 
have uniformly bounded variance. 

Theorem 1. Let the boundary condition hold with a smooth function cf). 


Let the Markov chain {XjjjgNo have the transition kernel P given by (18) and 
initial distribution IKq ~ 6^^, Xq G D\T, then the estimator (19) is convergent 
and unbiased. Moreover, there exits a constant C > 0, independent of Xq, 
such that 

Var[M(a;o, •)] < C'l 1^1 lioop^^)- 

Proof. To show that the sequence of successive approximations converges 
uniformly it is sufficient to prove the existence of a positive constant Ci < 1 


such that \\K^\ 


L°°{D\T) 


< Cl. For X ^ Dg we have 


/ / k{x, dy)k{y, dz) < max{l — pa{x)} < 1 

Jd\T Jd\T xedD 

and for x G D\{T U Zlg U T^) we may split the integral to obtain 


/_ /_ ^D\r(a^>d2/)A;(|/,dz) + /_ / P-p\Ti^,dy)k{y,dz) 

Id\t Jd\{tud^ut^) Jd\t Jd^ 


< 1 


*^ 0 ( 3 :, e) , ... ^o(a:,£) 

+ maxjl — Pa{x)\-^^^ -;rv;- < 1, 


2t: d-j^\rp[x) xedD 


‘^^d^\Tix) 


where we have used formula (10) specialized to the two-dimensional case and 


domain D\T. The convergence of the sequence of successive approximations 


yields convergence and unbiasedness of the estimator (19). 


Now let us write the estimator u{xo,uj) as the sum of local scores given by 
Sj{u;) := F{xj{u)) for all j < ((cu) and Sj(cj) := 0 for all j > ((cv). Obviously 
it holds that 

00 2 

Var[M(a:o,-)] < E^o[(5^Si(-)) • 

i=o 

By convergence of the successive approximations, the lifetime ( is P 2 ;g-a.s. 
hnite, implying that the series of local scores is absolutely convergent in 
square mean with respect to the probability space (O, Pa:(,). In particular 

we may write 


E, 


'XO 


00 

K5:' 

j=0 


= 2 




1=0 k=j 


1=0 
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By conditioning we obtain for j < k 

Exo[siOfc(-)] = < C] ■ ^xoU < 0 

and a straightforward calculation gives 




F{xj)F{xk)k{xo, dxi)...k{xk-i, dxk) 

J{D\T)>= 

= K^{FK^-^F){xo). 


Summation of these expectations yields 


OO OO OO OO 

EE K^{FK'^-^F) = ^ (f ^ K\f'^ ={I- K)-\F{I - K)-^F) 

j=0 k=j j=0 1=0 

and 

OO 

'^K^F^ = (I-K)-^F^. 

j=0 

Finally we arrive at 

Var|fi(x„,.)] < (I-A')-'(2F(I-A')-‘F-F2 )(i„) <C||F||J„,j 5 ^^). 


Notice that we have used the fact that by convergence of the successive 
approximation we may manipulate the Neumann series to obtain (Z—K)~^ = 
(X — + K) implying 

||(I-A)-'||i„,ii^^,<C2 = 2(l-Ci)-'. 

We have thus shown the assertion with C = C'2(2C'2 + !)• □ 

Let us conclude this section with an estimate for the mean square error 
of the partially reflecting RWOS estimator. 

Theorem 2. Let the boundary condition hold with a smooth function (j). 
For given £ > 0 there exist a step size h > 0 and a constant C > 0, such that 


E 



M 

lF{Xo,Um) 

m=l 
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Proof. The mean square error is equal to 


M 


M 


M 


E ( '^u%Xo,u}m) + (E[^M^(xo,a;m)] - u{xo) 


m=l 


Var[-u^(xo, ■)] 


m=l 

M 


m=l 


M 


+ e[E U%Xo,Um)] -u{X(i)) , 


m=l 


where the first term on the right-hand side is due to the Monte Carlo sampling 
error and the second term is due to the bias of the discretization. As the 
RWOS simulates the exit position exactly, the bias only comes from the finite 
difference discretization when hitting the boundary. It suffices to consider 
the case T = 0, as the variance of the estimator achieves its maximum in this 
case. Note that due to (12), e is required to achieve a local bias 0{h^). 


When the boundary is hit, the trajectory is absorbed with probability 
and the number of hits of the boundary follows a geometric distribution with 
this probability as parameter. The mean number of hits is thus given by 
1 -|- f. Consequently the global bias is of order 0{zhf). □ 


5. Approximation of electrode measurements 

As we have obtained a convergent partially reflecting RWOS estimator 
with uniformly bounded variance for the potential, we can immediately dehne 
an estimator for the electrode currents. Indeed we may write the potential as 
a sum u = Mo + CMi, where Uq and Ui solve auxiliary boundary value problems 
for the Laplace equation (|^ subject to the boundary conditions 

Mo = 0 on dT, zp ■ Vwolaz) + fuo\dD = 9 on dD, 

respectively. 


Ml = 1 on dT, zv ■ Vuilgci + fuilsD = 0 on dD. 
From the boundary condition ([^ one obtains 


U, 


= TIP I W ^ 


Mo(x) CMi(x) 




\^\ JEi Z 

and the conservation of charges condition ([^ yields 

N ^ .rr / N N 


dcT(x), I = 1,..., N, (21) 


Ui - Mo(x) 


1=1 


z 


N 

Mn) (EI 

1=1 P 


Ui{x) 


dcr(x) 


-1 
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Therefore we define for the integrals d(T(a;), i = 0,1 and I = 1, 

the estimators 


iW2 . Ml 

m2=l mi = l 




mi J ? 


( 22 ) 


using the so-called double randomization principle, cf. [ID]. That is, the 
potential is computed via the partially reflecting RWOS estimator and the 
boundary integrals in (21) are approximated via Monte Carlo sampling as 
well, using a uniform initial distribution Xq ~ lA{Ei), I = Hence 

A denotes an elementary element from the corresponding probability space 
and 

\E\ 

fi{xQ,u):= — -Ui{xo,u), i = 0,1. 

‘Z/ 

Obviously, to estimate the expectation, it would be sufficient to construct 
only one Markov chain for each realization of A. In practice, however, a 
splitting technique is usually used, where realizations of A and then for 
each of these realizations Mi independent Markov chains are constructed. 
For the optimal choice of Mi and M 2 we refer the reader to the book 


which yields the following estimator for the electrode currents: 


Convergence of the estimator (22) is an immediate consequence of Theorem 

(23) 


where each component Jf (Mi, M 2 , A,a;) is given by 


1 


U, 


_ j Mda(x)-^(ef,o(Mi,Mi,A,u;)+ef,i(Mi,Mi,A,u;)-c^ 


The constant c is approximated by the combined random estimator 


c" : = 


Ei=i ( Sb, i - ff,i)(Vi, M 2 , X, 01 
■EI,SAMuM2,X,oi) 


Remark 1. Note that in contrast to the idealized boundary condition the 
right-hand side of ^ presents some discontinuouities, so that the potential 
is merely Holder-continuous, cf. \4^ . In particular, the fourth order con¬ 
vergence with respect to h obtained in Theorem\^for the idealized boundary 
condition (flD will be reduced in the case of discrete electrode measurements. 


16 








6. Variance reduction 

In order to reduce the variance of the partially reflecting RWOS estimator, 
we propose a combined control variates conditional sampling technique. The 
basic idea of the control variates technique is to employ the known solution of 
a ‘nearby’ problem, see e.g. im. We shall exploit here the continuous depen¬ 
dence of the electrode currents on the conductivity. Let us thus consider the 
forward problem for the homogeneous medium with unit conductivity k = 1 
in D and let v denote the corresponding solution of the Laplace equation in 
D, subject to the boundary condition ([^. We may proceed as above and 
consider the inhomogenious integral equation 

v{x) = / v{y)k{x,dy) + F{x), x G D, (24) 

J D 

with the Radon measure 

{ (1 - Pa{.x))Ph{e^±e2){.x, dy), X e 

Pd(,x, dy), X e U T (25) 

P-^\T{x,dy), else, 

and inhomogeneity 


F(x) 


s(x) -|- hrQ^x), X E 
0 , else. 


As in the derivation of the partially reflecting RWOS estimator we dehne a 
Markov chain {XjjjgNo; this time, however, with state space (DU{c?}, Bd{D)) 
rather than {D\T U {5}, Bd{D\T)) and with transition kernel 


P{x,B):= I k{x,dy), xEDU{d}, B E Bd{D), (26) 

J B 

and initial distribution Xq ~ The key idea is that we may use one 
realization of this Markov chain to compute a realization of both, v^{xo, ■), 
as well as uf (xq, ■), f = 0,1. Let us dehne the practically realizable estimators 

1^1 , . 

r)f(xo,-) := ^(^n(xo) - h^(xo, •) +h-(xo,-)j, * = 0,1. 
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Let C denote the lifetime of then we obtain for i = 0,1, the condi¬ 

tional expectations 


^xo[v!{xo,-)\riT,)] 


^ f - v{Xr(T,)) + 0{e)), t{T^) < C 

^ [n(a;o), t{T^) > ( 


Now we set for i = 0,1, 


v!{xo, ■) := E^^[fil{xo, 


and taking the expectation yields for I = 1, ...,N and z = 0,1 

Ew(ijp[E^o(.)[r)-(a;o(-),-)]] = [ da{xo) + 0{e + h‘^). 

Jei z 

In particular the variance of fj^ixo, ■) is strictly smaller than the variance of 
fjl{xQ, ■) since 

Var[r7f(a;o, ■)] = Var[f}f(a;o, ■)] + ^xoi^^T^ivlixo, 


On top of that, we use n as a control variate. That is, either v is known 
explicitly, which is the case for certain geometries, cf. [I2l[36], or an approxi¬ 
mation of V via a hnite element or boundary element method is computed in 
a pre-computation step. In both cases we only need to simulate realizations 
of the random variable and then evaluate (27). In order to approxi¬ 

mate the electrode currents we proceed as above aimdehne for the integrals 
Ie dcr(a;), / = 1,..., iV, z = 0,1, the estimators 


^ Ma ^ Ml 

^2; A, cj) := ^ ^ ^ ^ hi (^o(-^m2); ^mi), * = 0,1, 

2 m2=l ^ ^ mi=l 


which yields a reduced variance estimator for the electrode currents if we 
substitute with - in the equations dehning (23). 


7. Generalizations 
7.1. Layered conductivities 

More realistic models in breast cancer modeling use layered conductivity 
models, see, e.g., |21]. Then the forward problem ([^, ([^ and ([^ is a diffrac¬ 
tion problem. Simulation of diffusion processes in piecewise constant media 
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has been studied in recent time, see, e.g. [2Sl ^ 1271 [2S]. The approach we 
adapt here was hrst introduced by Lejay and the hrst author in [25], therefore 
we content ourselves here with a brief description. For the sake of simplicity 
we assume that D is divided in two subdomains Di and D 2 := D\Di such 
that Di (Z D and T G Di. The interface S := dDi is assumed to be smooth 
and the conductivities in Di and D 2 will be denoted Ki and ^ 2 , respectively. 
We proceed similarly to the derivation of the partially reflecting RWOS es¬ 
timator, i.e., we use a hnite difference approximation in the interface layer 

S. 

h^A^u{xK + hei) = —hV^u{xK) + u{xk) — Rhuixx), (28) 

where we have assumed without loss of generality that I'i'XT.ixK)) = ei. The 
solution u of the diffraction problem is smooth in both subdomains Di\T and 
D 2 , continuous on S and satishes a transmission condition, i.e., the limit 

i^2{u{71y:{xk) + hei) - u{71t.{xk)) + i^i{u{71j^{xk) - hei) - u{71j^{xk))) 
h —^0 h 

vanishes. Let us introduce two parameters hi,/i 2 > 0, both of order 0{h), 
such that this transmission condition may be written in the form 

H2'^'1^u{71j:{xk)) = KiV"^W(7rs(a:if)) 0{K,2hl Kih\). 

In D\Di we obtain using the standard 5-point stencil 

H2hl/^’^^u{xK + ^ 261 ) = 0{K2hl) 

and in Di we have similarly 

Kih\^~’^^u{xK — hiCi) = 0{Kihl). 

Inserting those equations into ( |2^ yields 

K2u{xk) - K,2h2V’l^u{xK) - K2Rh2u{xK) = 0{K2h\) (29) 


and, respectively, 

Kiu{xk) — KihiV~^^u{xK) — K.iR-hiU{x k) = 0{Kih\). (30) 


Multiplying (29) by hi and (30) by h 2 and summing them up, one obtains 

Kih2 


u{xk) = 


K2hi 

K2hi + Kih2 


Rh^uixx) + 


K2hi Kih2 


R-hM^k) RrsixK), (31) 
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where t^^xk) = C>(k2^2^i + ^i^i^ 2 + £^) as hi, h2,£ —>■ 0 . As in the case of the 


partially reflecting RWOS estimator, expression (31) yields a probabilistic 
interpretation and thus a probabilistic estimator. By the strong Markov 
property of the partially reflecting Brownian motion one can couple this 
estimator accounting for the behavior at the interface S with the partially 
reflecting RWOS estimator. The resulting estimator may be analyzed in the 
same manner as described above for the partially reflecting RWOS estimator 
with constant background conductivity. 


7.1.1. Choice of the parameters 


We write equation (31) in the generic form 


u{xk) = pRh2u{xK) + (1 - p)R-hM^K) + rsixx), 


where p G (0,1). There are at least three natural choices of the parameters 
hi and h 2 , namely 


(i) h 

(ii) h 

(hi) h 


= K2hi 

= hi = 

_ h 

^ ~ 


= Kih 2 , then p = 1 — p = |, 
h 2 , thenp= A, 

h 2 = thenp = 


1 — p = 


K2+Kx ’ 

1 — P = 




y/i^+y/i^' 


Notice that the hrst choice is related to the kinetic scheme obtained in [28] . 
where the direction which was originally chosen uniformly in (0, 27r) is re¬ 
placed by a discrete random variable taking only 4 values. In (ii) the proba¬ 
bilities to go to one side of the interface correspond to those in [53] • Finally 
(hi) may be seen as a generalization of the one-dimensional scheme based on 
simulation of the skew Brownian motion in [25] . 

In our numerical examples, we will also have to deal with other types of 
boundary conditions and with multiple interfaces. Consequently, we chose 
method (ii) for all our tests to deal more easily with the constraints on the 
step h in order not to cross interfaces when replacing the motion. 


7.2. Uncentered walk on spheres (UWOS) 

When the physical domain is simple, one can compute the law of the 
exit point of the Brownian motion starting at any point of the domain ex- 
plicitely. This has been done for instance for rectangle domains in [TT] or for 
spherical domains in [SS] [31]. We restrict ourselves here to the case d = 2 
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and for a circle of radius R, centered at point (0, 0), the law of the exit po¬ 
sition of a Brownian motion starting at point (r cos(6'),r sin(6')) is given by 
(i?cos(a), i?sin(a)), where 

f R — T \ 

a \= 6 + 2 arctan (-tanfTrt/) ) 

\R + r ^ ’) 

and f/ is a uniform random variable in [0,1]. 

1.3. Random parameters 

In many practical situations the electrode currents depend on some ran¬ 
dom parameter /i, for instance due to random contact impedances, see e.g. 
[23]. In uncertainty quantihcation one is usually interested in calculating the 
expectation and the covariance of the random current measurements with 
respect to the law of this parameter. Computing these quantities in an ef- 
hcient manner is also crucial in the Bayesian modeling error approach, cf. 
|2Uj . where the statistical properties of modeling and discretization errors 
are estimated beforehand and subsequently used in the numerical solution 
of the inverse problem. As the underlying probability spaces are usually 
high-dimensional, uncertainty quantihcation suffers from the curse of dimen¬ 
sionality so that for many practical applications crude Monte Carlo sampling 
is still the method of choice. That is, one samples an ensemble of realizations 
of the random parameter and solves the deterministic boundary value prob¬ 
lem for each realization by a deterministic method such as a hnite element 
or boundary element method. However, the burden in terms of computation 
time of this procedure is likely to be prohibitive. 

In the framework presented here, this difficulty can be overcome natu¬ 
rally by using the double randomization principle, cf. mi, which yields the 
relations 

E^[JK0] 

Here uji,uj 2 and Ai,A 2 , respectively, are conditionally independent trajecto¬ 
ries constructed for a hxed realization of /i. 

8. Numerical tests 

Our numerical tests were performed using a circular model, i.e., D is the 
planar unit circle, see Figure Such a geometry may serve as an appropriate 
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Figure 3: The benchmark setting modeling a breast geometry with 8 electrodes and a 
layered conductivity. The electrodes are numbered clockwise. 


model for certain mammography systems, cf., e.g., mm- We chose a circular 
inclusion of radius r centered at the origin. In this case the constant on 
the boundary of the inclusion must be equal to zero for symmetry reasons. 
We used 8 electrodes, each of width 0.1. An alternating voltage pattern 
was imposed, i.e., Uj = (—l)b j = 1,...,8. For the computation of the 
electrode currents, we used the double randomization principle, that is, the 
starting point of each trajectory was picked uniformly at random on one of 
the electrodes. 

The numerical scheme was implemented in Fortran and parallelized via 
OpenMP. The test cases were run on a workstation with 4 AMD Opteron 8 
core CPUs. Pseudo random numbers were generated with an implementation 
of L’Ecuyers’s parallel MRG32k3a random number generator. The reference 
solution was computed using hnite element routines from the EIDORS pack¬ 
age, cf. [45] . 


8.1. Unit background conductivity 

8.1.1. Idealized measurement model 

To verify the theoretical result of Theorem 2 numerically, let us hrst 
study the idealized measurement model assuming that measurements can be 
taken on the whole boundary. In our experiment, we considered the Robin 
boundary condition 

zVuldD + u\aD = cos(46'). 


where 6 denotes the polar angle. In order to analyse the global bias of the 
partially reflecting RWOS estimator, we computed the bias of the 

practically realizable estimator (20) at a single point x = (0.99361,0.11286) 
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Figure 4: Idealized measurement model: The logarithm of the bias is plotted against the 
logarithm of the stepsize; solid lines show the corresponding least-squares fits. 

for two different values of the contact impedance z, namely 

Zi := 0.5, Z2 := 0.1 

and 5 different values of h chosen equidistantly from the interval [0.08,0.2]. 
The reference values for a centered circular inclusion of radius r = 0.3 were 
computed using a very hne discretization by linear hnite elements, cf. |1^; 
we found u{x) ~ 0.299 and u{x) ~ 0.642, respectively. In Figure]^ we plot 
for each contact impedance the bias corresponding to the approximations 
obtained by the new estimator using 10® simulations and the different step- 
sizes in a logarithmic scale together with the corresponding least-square hts. 
We obtained an estimated order of convergence (EOC) of 2.05 and 2.07 for 
Zi and Z 2 , respectively. Moreover, as one would expect from the proof of 
Theorem 2, the bias increased, when ^ increased. 

8.1.2. Discrete electrode measurements 

Bias estimation. In this experiment, we computed for the direct method (i.e. 
the method without variance reduction) the bias B^., i = 1,2, of the electrode 
current through the electrode centered at (1,0) for 5 different values of h 
chosen equidistantly from the interval [0.04,0.2]. 

Again, the reference values for a centered circular inclusion of radius r = 
0.3 were computed using a very hne discretization by linear hnite elements. 
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Figure 5: Discrete electrode measurements: Finite element discretization and the reference 
solution within the domain D for an alternating voltage pattern. 

see Figure In Figure we plot for each contact impedance the bias 
corresponding to the approximations obtained by the new estimator using 
10 ® simulations and the different stepsizes in a logarithmic scale together 
with the corresponding least-square fits. We obtained an EOC of 1.76 and 
1.62 for zi and ^ 2 , respectively. As one would expect from Remark the 
EOC is reduced compared to the idealized measurement model. 

Variance reduction. Next, we investigated the efficiency of the variance re¬ 
duction method proposed in Section]^ The control variate, i.e., the electrode 
current corresponding to the problem without inclusion was precomputed us¬ 
ing a finite element discretization. Table [T] shows the experimental results for 
the test problem described in the previous paragraph with contact impedance 
Z 2 and an inclusion centered at the origin with different radii. 


r 

Tref 

•^3 

<^3 

’^3 

<^3 

0.9 

0.976 

0.974 

0.445 

0.974 

0.161 

0.8 

0.902 

0.896 

0.517 

0.904 

0.098 

0.7 

0.874 

0.868 

0.612 

0.876 

0.054 

0.5 

0.864 

0.872 

0.681 

0.865 

0.015 

0.3 

0.862 

0.871 

0.751 

0.862 

0.001 


Table 1: Variance reduction: Approximations of the reference electrode current via 
the direct method J| and via the method with variance reduction J| together with the 
corresponding standard deviations. 
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Figure 6: Discrete electrode measurements: The logarithm of the bias is plotted against 
the logarithm of the stepsize; solid lines show the corresponding least-squares hts. 

The reference values for the different radii of the inclusion were com¬ 
puted using a very hne discretization by linear hnite elements; All simulations 
were performed using 10® Monte Carlo simulations with stepsize h = 0.004 
and e = 10“®. We computed approximations of the electrode current denoted 
J| and its standard deviation df for the direct method, respectively and 
(j| for the method with variance reduction. 

As one would expect, the standard deviation increases when r decreases 
for the direct method. For the method with variance reduction, the standard 
deviation decreases when r decreases because the problem gets ‘closer’ to 
the one without inclusion which is used as control variate. This means in 
particular that the efficiency of the method with variance reduction increases 
as the size of the inclusion decreases which is particular interesting with 
regard to the inverse problem. 

Subsequently, we investigated the efficiency of our approximation meth¬ 
ods, as well as different variants thereof, based on the quantity C given by 
the variance multiplied by the computational time which is a standard cri- 
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terion for Monte Carlo methods. Note that for the variants of the method 
with variance reduction, which shall be described below, we did not include 
the computational time required to solve the problem without inclusion as 
this is a precomputation which can be done once and for all. The results 
for the direct method, the method with variance reduction based on FEM, 
the RWOS and the UWOS, respectively, are shown in Table The latter 
methods use the RWOS, respectively the UWOS instead of the FEM to com¬ 
pute the solution on the inclusion. For these methods, the superscript in the 
notation denotes the number of sample paths used for this computation. We 
see that the method with variance reduction based on FEM is more efficient 
than the other variants and that its efficiency increases as r increases. 

8.2. Discrete electrode measurements and layered background conductivity 

Now let us turn to the more realistic case of a layered background con¬ 
ductivity. In this experiment, we considered a domain D which is separated 
into two areas with diffusion coefficient Ki for all points (x, y) such that 
y/x^ -f- > R and with diffusion coefficient K 2 elsewhere outside the inclu¬ 

sion, see Figure]^ We chose ki = 1.5, ^2 = 1 and R = 0.9. 

8.2.1. Bias estimation 

We computed the bias i = 1, 2, of the of the electrode current through 
the electrode R3 centered at (1,0) for 5 different values of h chosen equidis- 
tantly in the interval [0.02, 0.1]. As in the previous experiments, the reference 
values for a centered circular inclusion of radius r = 0.3 were computed us¬ 
ing a very hne discretization by linear hnite elements. In Figure we plot 
the bias corresponding to the approximations obtained by the direct method 
using 10® simulations in a logarithmic scale together with the corresponding 
least-square hts. 


r 

Cj^ir 

Cfe 

'^RW 

'^RW 

^(10) 

'^RW 

'^UW 

'^UW 

'^uw 

0.9 

0.74 

0.78 

16.7 

16.2 

19.5 

5.5 

4.1 

3.7 

0.8 

2.0 

0.21 

14.2 

12.5 

12.3 

6.1 

4.0 

2.5 

0.7 

3.4 

0.07 

13.0 

9.9 

8.4 

6.0 

3.9 

1.9 

0.5 

6.6 

3.0E-3 

9.1 

7.2 

4.8 

5.9 

3.5 

1.3 

0.3 

10.1 

4.5E-5 

7.2 

4.8 

2.9 

5.5 

2.9 

0.9 


Table 2: Efficiency of the methods (unit background): The direct method, the method 
with variance reduction based on FEM, the RWOS and the UWOS, respectively. 
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Figure 7: Discrete electrode measurements and layered background conductivity: The 
logarithm of the bias is plotted against the logarithm of the stepsize; solid lines show the 
corresponding least-squares fits. 

We obtained an EOC of 1.80 and 1.86 for zi and Z 2 , respectively. In 
analogy to the observation in Remark [T| the discontinuity of the conductivity 
leads to a merely Holder-continuous so that the EOC is slightly smaller than 
two. 

8.S.S. Variance reduction 

Next, as in the previous experiment, we investigated the efficiency of the 
variance reduction method proposed in Section 6. The control variate, i.e., 
the electrode current corresponding to the problem without inclusion was 
precomputed using a hnite element discretization. As in the previous exper¬ 
iment, we investigated the efficiency of our approximation methods as well 
as different variants thereof, based on the quantity C given by the variance 
multiplied by the computational time. The results for the direct method, the 
method with variance reduction based on FEM, the RWOS and the UWOS, 
respectively, are shown in Table The latter methods use the RWOS, re¬ 
spectively the UWOS instead of the FEM to compute the solution on the 
inclusion. For these methods, the superscript in the notation denotes the 
number of sample paths used for this computation. Note the increase of the 
values of C compared to the previous method for the direct method. This 
can be explained by an increase of the computational times due to the time 
spent by the trajectory near the interface. Moreover, we see that the variance 
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r 

C]Jir 

Cfe 

/^W 

^RW 

^RW 

'^RW 

^UW 

^UW 

^(10) 

'^uw 

0.8 

4 

0.27 

49 

45 

43 

44 

40 

38 

0.7 

8 

0.09 

45 

39 

34 

41 

34 

31 

0.6 

12 

0.03 

42 

33 

27 

38 

29 

24 

0.5 

17 

8E-3 

31 

25 

20 

30 

23 

17 

0.3 

28 

9E-5 

29 

22 

15 

28 

20 

12 


Table 3: Efficiency of the methods (layered background): The direct method, the method 
with variance reduction based on FEM, the RWOS and the UWOS, respectively. 

reduction method based on FE is still highly efficient and clearly superior to 
the direct method. 

8.3. Random background conductivity 

In many practical situations, the background is not perfectly known and 
is therefore modelled as a random medium. We consider a layered model 
where the diffusion coefficient in each region is given by uniformly distributed 
random variable. We assume that ki is uniformly distributed in the interval 
[1.3,1.7] and /t 2 is uniformly distributed in the interval [0.8,1.2]. We also 
assume that the radius i? is a random variable, uniformly distributed in 
[0.89, 0.91]. As described in Section]^ the computation of the mean value of 
the electrode current is computed using the double randomization principle. 
A realization of the medium is picked according to its distribution and a 
starting point for the trajectory is picked according to a uniform distribution 
on the electrode; then one trajectory is simulated and the corresponding 
score is computed. This procedure is repeated a number of times and the 
resulting scores are hnally averaged. Concerning the variance reduction, it 
is not obvious how to use a FEM approximation of the control variate in a 
straight forward way without remeshing for each realization of the random 
medium. Therefore, we restrict ourselves here to the control variate based on 
the continuation of the walk via the UWOS method with 10 trajectories. We 
note that the approximation of the mean value of the solution to the partial 
differential equation with a stochastic coefficient is computed without extra 
cost and with a variance of the same order as the one for the deterministic 
problem. This is a huge advantage compared to deterministic methods, where 
the solution of an elliptic boundary value problem is required for each draw 
corresponding to a realization of the random medium. In our experiments, 
the probabilistic methods were up to 100 times faster than the Monte Carlo 
















r 



<^3 

COir 

^3 

^3 

'^uw 

0.7 

0.902 

0.896 

0.654 

7 

0.902 

0.225 

28 

0.6 

0.874 

0.874 

0.703 

12 

0.872 

0.207 

22 

0.5 

0.864 

0.860 

0.747 

17 

0.865 

0.153 

19 

0.3 

0.862 

0.860 

0.830 

28 

0.860 

0.102 

12 

0.2 

0.862 

0.860 

0.899 

36 

0.863 

0.057 

9 


Table 4: Variance reduction and efficiency of the methods (random background): Ap¬ 
proximations of the reference electrode current via the direct method J| and via the 
method with variance reduction Jf together with the corresponding standard deviations 
and the quantity C. 

sampling based on FEM solution of the forward problem. Moreover, the 
variance reduction is still efficient to some extent which is an important 
aspect with regard to the inverse problem of detecting small anomalies. 

9. Summary and future work 

The complete electrode model is the most realistic model to approximate 
real electrode measurements in electrical impedance tomography. We have 
given a probabilistic interpretation of the corresponding electrode currents 
taking into account both the mixed boundary condition and the possibly 
discontinuous diffusion coefficient. Then, we have proposed a Monte Carlo 
method based an a novel partially reflecting RWOS estimator to compute 
these currents in an embarrassingly parallel manner. This method involves 
the RWOS algorithm inside subdomains where the diffusion coefficient is con¬ 
stant and replacement techniques motivated by hnite difference discretization 
to deal with mixed boundary conditions as well as transmission conditions. 
The global bias of the corresponding algorithm is analyzed both theoretically 
and experimentally. Moreover the variance of the new estimator is stud¬ 
ied and subsequently considerably reduced via a control variate conditional 
sampling technique. Indeed, it is this variance reduction which makes the 
proposed method such an interesting alternative to standard deterministic 
methods, even for two-dimensional problems. 

In future work, we intend to use the new Monte Carlo method in the 
framework of Bayesian statistical inverse problems. We expect that the in¬ 
herent parallelism of our method, combined with the wide availability of 
multi- and many-core computing hardware, will enable the efficient treat- 
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ment of three-dimensional problems which would be prohibitively expensive 
in terms of computation time with a (non-parallelized) deterministic forward 
solver. 
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